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Gaussian random fields pervade all areas of science. However, it is often the departures from 
Gaussianity that carry the crucial signature of the nonlinear mechanisms at the heart of diverse 
phenomena, ranging from structure formation in condensed matter and cosmology to biomedical 
imaging. The standard test of non-Gaussianity is to measure higher order correlation functions. In 
the present work, we take a different route. We show how geometric and topological properties of 
Gaussian fields, such as the statistics of extrema, are modified by the presence of a non-Gaussian 
perturbation. The resulting discrepancies give an independent way to detect and quantify non- 
Gaussianities. In our treatment, we consider both local and nonlocal mechanisms that generate 
non-Gaussian fields, both statically and dynamically through nonlinear diffusion. 



Random fields pervade all areas of science. A disparate 
class of phenomena, ranging from the cosmic background 
radiation p] and surface roughness [2] to medical images 
of brain activity [5| and optical speckle patterns [3] , pro- 
duce data that can be regarded as random fields. The 
statistics of geometrical features of these fields, such as 
the density of extrema of various types, can be used to 
characterize them [21 E] . When the fields can be approx- 
imated as Gaussian fields, the physical meaning of these 
statistical properties is generally well understood [S] : 
the statistics of extrema reflects the amount of field fluc- 
tuations at short distances. 

Though analytical investigations are often restricted 
to Gaussian fields, phenomena described by nonlinear 
laws (such as the dynamics of inflation that produced the 
cosmic background radiation) produce non-Gaussian sig- 
nals. Quite often, the observable signal is averaged over 
a large scale, producing approximately Gaussian statis- 
tics on account of the central limit theorem, and masking 
the nonlinearity. Nevertheless, the surviving tiny depar- 
tures from Gaussianity can carry the crucial signature of 
the nonlinear microscopic mechanisms at the heart of the 
phenomena. As an illustration, consider a low resolution 
measurement of the spatial magnetization of a material 
well above the critical temperature. The magnetization 
fluctuates like a Gaussian random variable - each region 
contains many domains oriented up or down in arbitrary 
proportion. However, a small non-Gaussian contribution 
remains, because there is a maximum possible magne- 
tization per unit area which can be traced all the way 
down to the quantization of the spin of the electrons and 
hence the probability distribution cannot exhibit Gaus- 
sian tails. 

In order to unveil such elusive effects, one needs an 
indicator that is sensitive to both short distances and 
small signals. 

The most common tool used to probe the statistics 
of a random field is to measure its correlation func- 
tions. For example, the statistical properties of a ran- 



dom scalar field, h(x,y), with Gaussian statistics, are 
entirely determined by its two-point correlation function 
(h(x, y)h(x', y')), and its higher-order correlation func- 
tions can be written simply as the sum of products of 
two-point correlation functions. The nonfactorizability 
of these higher-order correlation functions is one of the 
standard indicators of non-Gaussian statistics. 



Here we focus on a more geometric approach: view the 
scalar field as the height of a surface and study its ran- 
dom topography to infer the statistical properties of the 
signal (see inset of Fig. [I]). The densities of peaks and 
valleys, or of topological defects in the curvature lines 
known as umbilics (see Fig. [TJ , are sensitive indicators of 
how jagged the height field is at short distances. As we 
shall see, they provide an independent pipeline to detect 
non-Gaussianities, distinct from multiple-point correla- 
tion functions. This geometric approach has been applied 
successfully to track the power spectrum of a Gaussian 
field and it has been the subject of extensive theoretical 
and experimental studies [oT-fTU]. 

In this paper, we introduce the key physical con- 
cepts and mathematical techniques necessary to study 
the stochastic geometry of signals that can be described 
as a Gaussian random field plus a perturbation that 
we wish to track. We first show how to treat non- 
Gaussianities within a local approximation and calculate 
how the statistics of extrema change when a nonlinear 
transformation Fml{Hq) is applied locally to a Gaus- 
sian field Hq- Then we consider the case of fields that 
cannot be probed directly, by calculating the statistics 
of umbilical points, which are topological defects of the 
lines of principal curvature Finally, we turn to the 
class of nonlinear diffusion equations [2 and go beyond 
the local approximation, by considering the effects of spa- 
tial gradients that couple values of the field at different 
locations. As an illustration, we solve explicitly for the 
nonlocal non-Gaussianities generated dynamically by the 
deterministic KPZ equation which models surface growth 
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Figure 1: The principal curvature of a surface. The major and 
minor axes of the small ellipses represent the direction and rel- 
ative size of maximum and minimum curvature at the centre, 
equivalent to the direction of maximum and minimum po- 
larization for an optical field. The curvature lines are always 
tangent to the direction of maximal curvature. At some points 
the curvature is the same along all directions (the equivalent 
polarization is circular); these are called umbilical points, of 
which there are three types, all shown in this image: on the 
top left is a lemon, on the top right a star and on the bottom 
right a monstar. The circles demonstrate their topological in- 
dices: — 1/2 for a star, +1/2 for the other two. The lemon 
has one (locally) straight curvature line terminating at it (in- 
dicated with a thick line), the other two have three. The inset 
shows a computer-generated Gaussian surface with periodic 
boundary conditions, a small square of which served as the 
source of this picture. 

I. CRITICAL POINTS 

To gain some insights into the physical mechanisms 
that generate non-Gaussianities, consider first how an 
isotropic Gaussian field Hc{r) arises from the random 
superposition of waves (or equivalently, Fourier modes) 

H G (r)=J2 A (k)cos(k-f+^), (1) 
% 

with an amplitude spectrum A(fc)|21j that depends only 
on the magnitude of the wave vectors, k = |fc|. The 
phases 0j are uncorrelated random variables uniformly 
distributed in the range [0, 2-k\. The statistical proper- 
ties of the Gaussian field Hair) are entirely encoded by 
the function A(k), or equivalently, its moments Km — 
J2 K k 2n iA(k) 2 . 

The most basic difference between Gaussian and non- 
Gaussian variables is that Gaussian ones are always sym- 



metric about their mean. As a consequence, irrespective 
of its power spectrum, a Gaussian field has equal den- 
sities of maxima and minima. Hence, a nonvanishing 
imbalance An between these two types of extrema serves 
as a probe to detect and quantify the non- Gaussian com- 
ponent of a signal, provided that it can be measured di- 
rectly. 

For example, consider the primordial curvature pertur- 
bation field, <&, a nearly Gaussian field of central inter- 
est to modern cosmological studies pp. Within a local 
approximation, the primordial field is obtained from a 
Gaussian field $g via a nonlinear relation $ = &q + 
fni^a +9m®G- Determining the parameters /„/ and g n i 
is one of the central tasks in the study of cosmological 
non-Gaussianities. As we shall see, the quadratic coeffi- 
cient can be determined from the imbalance An between 
maxima and minima of <fr. 

The imbalance can be derived in the more general con- 
text of a non-Gaussian field h that is obtained from a 
Gaussian field Hq via any nonlinear deformation h = 
Fnl(Hg)- If Fnl is a monotonic function, the maxima 
and minima do not change - only a nonmonotonic be- 
havior of Fml can alter this balance. 

The critical points of h are given by V/i = 
F' nl {Hq)V Hq = 0, where the dash indicates the deriva- 
tive of Fnl(Hg) with respect to Hq. This condition 
shows that h and Hq have the same critical points. Note 
however that, if F' NL (HQ(fo)) < at a critical point rb, 
then a maximum (minimum) at Hq^q) will be turned 
into a minimum (maximum) at h{fo). The number of 
saddle points does not change because of topological con- 
straints (it is equal to the invariant sum of maxima and 
minima). 

If the transformation has a bias towards converting 
minima into maxima, then h will have more maxima than 
minima; for example, h = Hq +eH g reverses its slope at 
sufficiently negative values of Hq, which are most likely 
to be minima. Following this logic, the first step toward 
calculating the imbalance between maxima and minima 
is to determine the probability g(z) that flg(ro) = z 
for a minimum ro of Hq. The symmetry properties of 
Hq imply that the analogous probability distribution for 
maxima is g{— z). 

The fraction of minima of Hq that become maxima of 
h is obtained by integrating g(z) over the range of z for 
which F' NL {z) < 0. Likewise, the fraction of maxima of 
Hq that are turned into minima is given by the integral 
of g{—z) over the same range. The overall imbalance 
in the densities of the maxima and minima of h can be 
readily obtained by adding these opposite contributions. 
The result reads 

A^max ^min 
n = 

^max T" n m i n 

= dz (g(z) - g(-z)). 

Jz:F' NL (z)<0 

For two dimensions, the exact analytical expression for 
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Figure 2: The relative difference An between the densities 
of maxima and minima of h = Hq + £Hq, where Hq is a 
Gaussian field with A = 3/4, as a function of e. The data 
points are results from computer-generated fields, the solid 
curve is the theoretical result (Eq. Q). The inset shows the 
corresponding distribution of minima g(z) (which forms the 
basis of our theoretical result), both the theoretical curve and 
a histogram of data gathered from computer-generated fields. 



g{z) is explicitly derived in Appendix A - it depends only 
on the moments Kq, K2 and K±. Rescaling h(f) — > h(ar) 
does not affect the function g(z), since it increases the 
density of maxima with any value of Hq by the same 
proportion. Hence, only Kq — (H G ) (which sets the scale 
of the distribution) and the dimensionless parameter A = 

A 2 
A',, A', 



can enter in the expression for g(z) (see plot in the 
inset of Fig. [2). 

As an illustration, apply Eq. ^ to the perturbed 
Gaussian field h = Hq + sH G . Figure [2] shows our theo- 
retical formula as a continuous line, validated by numer- 
ical data (dots) obtained from computer-generated ran- 
dom surfaces with amplitude spectrum A(k) ~ 6{kr>~k), 
having A = |. 

The imbalance between maxima and minima is partic- 
ularly useful to track large deviations from Gaussianity. 
This can be seen explicitly for the quadratic perturba- 
tion considered above. Note that, since Hq itself has an 
equal number of maxima and minima, the imbalance An 
in h is only created when one of these critical points is 
inverted. Whether Hq has a high likelihood of having a 
negative F' nl (Hq) = 1 + 2eHq is controlled by e^J {Hq). 
The probability of \2eHq\ exceeding 1 is exponentially 
small. Thus, An rises roughly as e~ c ^ e ' i? o)) 5 where 
c s» |. Note that our approach in deriving Eq. ^ and 
g(z) is nonperturbative and hence capable of describing 
large deviations from Gaussianity (for the specific type 
of deviation considered here), as demonstrated by the 
agreement between our formula and numerics over the 
entire e-range probed in Fig. [2j 



II. UMBILICAL POINTS 



The near-Gaussian fields under investigation are not 
always directly accessible experimentally. For example, 
the mass distribution along the line of sight responsible 
for weak gravitational lensing is mostly composed of dark 
matter and hence it cannot be detected directly [T2]. If 
the projected gravitational potential over a flat patch of 
the sky is taken to be the height of a 2D surface [13], the 
measurable shear field, is given by the lines of principal 
curvature [TT], as shown in Fig.JT] At some special points 
called umbilics, the curvature is equal in all directions, 
so the shear field cannot be defined and it must vanish. 
More precisely, a point r= ona surface with height 

function h(f) is an umbilic if the second derivatives sat- 
isfy the two conditions h xx {r) = h yy (r) and h xy {r) = 0. 
The ratio between different types of umbilical points 
(which is a universal number for an isotropic Gaussian 
field) serves as an indicator of non-Gaussianities in lieu 
of the extrema which cannot be detected. A similar rea- 
soning can be applied to study polarization singularities 
in the cosmic microwave background 14J and topological 
defects in a nematic or superfluid near criticality [151 116j . 

Inspection of Fig. [I] reveals that there are three types 
of umbilics: lemons, monstars and stars. Note that these 
umbilics are topological defects in the curvature- line field. 
The topological index of any umbilic is equal to plus 
(minus) 1/2, if the curvature-line field rotates clockwise 
(counter-clockwise) by an angle 7r, along any closed path 
encircling that umbilic only in the clockwise direction. 

A star has three curvature lines terminating at it and 
a topological index of — |. A lemon has only one line 
and index +5. A monstar has index +5, like a lemon, 
but three lines terminating at it, like a star. A striking 
feature of isotropic Gaussian fields is that the monstar 
fraction, the relative density of monstars with respect 
to all umbilics, equals olm = \ — = 0.053; this is 
a universal number independent of the power spectrum 
[U [TU] . Any deviation from this special value is therefore 
a sure sign of non-Gaussian effects. 

We will again consider a height field that is a nonlinear 
function of a Gaussian, h = F^^(Hq). In this general 
case, we were not able to find an exact result as we were 
for the extrema. We shall therefore assume that the non- 
linear contribution is small, such that h = HQ + ef(HQ), 
with / a nonlinear perturbation and e « 1 a small pa- 
rameter controlling the size of the nonlinearity. We can 
express cxm in terms of the joint probability distribution 
p of the second and third derivatives of h. To calculate 
p, we use the property that a probability distribution is 
determined by its moments of all orders (if the distribu- 
tion is well-behaved). This is most easily accomplished 
by calculating the generating function of the distribution, 
X - its logarithm can be directly expressed in terms of 
the cumulants C n 
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Table I: All nonzero cumulants 



C 2 (/i zz , h z * z » ) 


a(l + 2(/ (H))) 


Clihzzz, h z *z*z*) 


r(l + 2(/'(fl)» 


C2(hzzz* , hzz-z* ) 
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Cz(hzz, hzzz* , h Z 'z*z*) & conj. 


-3a 2 (/"(//)) 
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k conj. -6a 3 </"'(#)) 



logx(Ai,...,A„) 



E 



Cm(^jl | • • • ) Oin ) Aji 



(3) 

where £1, • • • , £n are the stochastic variables given by the 
spatial derivatives of h. The probability distribution can 
then be obtained from the generating function by taking 
the inverse Fourier transform with respect to Ai, . . . , A„. 
The cumulants can be written in terms of expectation 
values, e.g. C 2 (£i,6) = - In tnis context, 

these expectation values are called moments, which are 
not to be confused with the moments Km defined before. 

For Gaussian variables, only the second-order cumu- 
lants are nonzero. This gives rise to a generating function 
of the form 



(4) 



One can easily check that the inverse Fourier transform 
of x hi Eq. Q precisely yields the probability distribu- 
tion for a set of correlated Gaussian variables (assuming 
(d) = 0), see Eq. (Al). More generally, by determining 
all the cumulants, one can construct the generating func- 
tion and from that the probability distribution. We shall 
derive the monstar fraction up to first order in the pertur- 
bation ef(Hc) only; consistently we need to determine 
all the cumulants up to first order only. 

The monstar fraction (even at order s) could in princi- 
ple depend on / in a complicated way, if / is an arbitrary 
nonlinear function. In fact, quadratic terms in the func- 
tion / produce degree 3 cumulants in the distribution 
function of the field h (i.e., skewness), cubic terms pro- 
duce kurtosis (degree 4 cumulants) and in general degree 
n terms in / produce degree n + 1 cumulants. How- 
ever, the monstar fraction can be determined from just 
the distribution of a few derivatives of h, whose cumu- 
lants vanish beyond the fourth order due to symmetry, 
as shown in the Supporting Information. Consequently, 
the final result for the monstar fraction depends only on 

a single parameter, (f"(H G )) = f"(u)e^> 7 §=, 
where the primes indicate derivatives with respect to Hq. 

The calculation can be briefly summarized as follows. 
The monstar fraction is related to the distribution func- 
tion of some of the second and third derivatives of h, 
which we write in terms of complex coordinates z = x+iy 
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Figure 3: The fraction of monstars (see inset) olm of h = 
Ha + eHq, where Hg is a Gaussian field with fx = 16/27, as 
a function of e. The data points are results from computer- 
generated fields, the solid line is the theoretical first-order 
result (Eq. |5|). At e = we retrieve the universal fraction 
olm — 1/2 — l/v5 = 0.053, valid for any isotropic Gaussian 
field. 



and z*\ h zz , h zzz , and h zzz *. The definition of an um- 
bilic point becomes h zz = 0, where h zz is now complex. 
All the cumulants of these variables and their conjugates 
may now be calculated (up to order e). The complex 
coordinates allow for optimal usage of rotational and 
translational symmetry. Only a few of the cumulants 
are nonzero and these are evaluated in Table |TJ With 
the aid of these cumulants, the generating function can 
be constructed to first order using Eq. ^. Taking the 
Fourier transform leads to the probability distribution 
p(h zz , h ZZZl h zzz * ); the explicit form is rather long, but 
it is basically a Gaussian perturbed by cubic and quartic 
terms in h and its derivatives [20j . To obtain cvm, we 
then have to set h zz — and integrate over h zzz and 
h zzz * (taking care to include the appropriate Jacobian 
factor). Integrating over all of C 2 gives the total density 
of umbilical points, while the density of monstars is ob- 
tained by integrating over a specific range, which can be 
found in Appendix B. The monstar fraction is then the 
ratio of these two densities. The resulting deviation from 
am = 0.053 is 



Aa M = 0A29n(f'"(H G ))e, 



(5) 



where /i = K^/Kq. When applied to the local modal of 
the primordial field $ described before, Aaj\f in Eq. ^ 
depends only on the cubic coefficient g n i and not on 

Hence, the leading order perturbation that alters the 
monstar fraction is /(Hg) = H G , for which (f" (Hg)) = 
6. Note that this perturbation, like any odd and/or 
monotonic function of Hq, does not have an effect on 
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Figure 4: The imbalance between maxima and minima An, 
as a function of time, for an initially Gaussian field evolving 
according to the deterministic KPZ equation (Eq. (|6|). At 
t = 0, the Gaussian field was taken to have a Gaussian power 
spectrum (A(k) 2 ~ exp(— fc 2 /2fco)); see inset on the left. As 
time evolves, the surface becomes smoother (see inset on the 
right), decreasing the densities of maxima and minima, but 
also creating an imbalance between the two. The data points 
stem from simulations, for which A/4i/ = 0.1 was used. The 
solid curve is the theoretical result. 

the density of maxima and minima. Figure [3] shows 
ctM, as determined by Eq. ^ (continuous line), to- 
gether with data from computer simulations (symbols): 
the agreement between theory and numerics is very good 
in the linear regime. The spectrum used was again 
A(k) ~ 9{kn — fc), for which fi = Note that the mon- 
star fraction is very sensitive to a small non-Gaussianity, 
with a 20% change when e is just 0.01. For larger values 
of e, nonlinear effects become important and prevent ajj 
from becoming negative. In this regime, our approximate 
result must break down. 



III. NONLOCAL MODEL AND EVOLUTION 
EQUATIONS 

In section ||J we presented an exact expression for the 
imbalance between maxima and minima for a local per- 
turbation of a Gaussian random field. This simple class of 
models may describe the local evolution of a system that 
starts out with a Gaussian distribution, as for instance 
the growth of a population of cells that are initially dis- 
tributed on a dish and then divide without any significant 
migration from one region to another. 

However, many dynamical systems evolve in a non- 
local, nonlinear way. Non-Gaussianities are generated 
dynamically from the nonlinear equations of motion 
that the field h(r, t) obeys, even if the initial condition 



h(r,0) — Hc{r) is Gaussian. Unlike the case of local 
evolution, the imbalance between maxima and minima 
will now exhibit a power law increase, as the nonlinear 
perturbation grows. 

A broad class of nonlinear diffusion equations describes 
the necessary mixing between regions. Examples include 
several models of structure formation in both condensed 
matter [T7] and cosmology pQ , the Cahn-Hilliard equation 
for the development of order after a phase transition [TS] 
and simplified models of surface growth [5] . We will focus 
on the last of these, the deterministic KPZ equation [2] 
which models the height evolution of a substrate as atoms 
accumulate on it: 



U ^=vV*h{r,t) + ^hmf- (6) 

This equation is arrived at in the following way: to first 
order, the surface simply grows at a constant rate. This 
constant rate does not appear in the equation, however, 
since it is simply subtracted out; the two terms on the 
right-hand side are the additional effects. The first term 
describes the diffusion of particles along the surface, and 
the second nonlinear term describes approximately how 
the growth-rate varies with the local slope. The surface 
is assumed to grow at a constant rate perpendicular to 
itself, but since the height is measured vertically, h de- 
pends on the slope: this gives rise to the term quadratic 
in V/i. 

Our approach can be applied to other nonlinear dif- 
fusion problems well beyond surface-growth dynamics, 
when a different choice for the quadratic term is made. 
For example, if the term quadratic in the gradient in 
Eq. ^ is substituted by a term quadratic in the field, 
— h 2 , one obtains the Fisher equation which describes the 
growth and saturation of a population. A third possi- 
bility is to consider a mixed term hVh which gives the 
Burgers' equation governing shock dynamics and traffic 
flow. In all of these cases, we can study the time evolu- 
tion of an initially Gaussian h(f) profile. Upon setting 
the coefficient of the nonlinear term A equal to zero, we 
always retrieve the heat equation, which preserves the 
Gaussianity of h for all later times. On the other hand, 
for A 7^ 0, h becomes non-Gaussian. 

For concreteness, we discuss how an imbalance be- 
tween maxima and minima is generated by nonlocal 
non-Gaussianities in the context of the KPZ equation 
(Eq. ([6])). The nonlinear term breaks the symmetry be- 
tween positive and negative values of h which is a neces- 
sary condition to generate an imbalance. Note however 
that, in the case of a local evolution, this imbalance was 
exponentially small because the local evolution cannot 
create new extrema - it can only convert a maximum 
into a minimum whenever h happens to have a sufficiently 
large fluctuation. It is the presence of the diffusion term 
that is able to create new maxima and minima, even 
though, on its own, it would not be able to generate any 
imbalance, because of the symmetry h — > —h. The two 
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Table II: All second and third order cumulants 



<r=(\h z \ 2 ) 

« = <IM 2 > 



7 = 
5 = 



\hi 

u3 



\h zz *) 

{hi..) 
\h zz \ 2 h zz *) 



terms on the right-hand side of Eq. (|6| conspire together 
to change the number of maxima and minima asymmet- 
rically. 

Since the imbalance between maxima and minima will 
now have a contribution pcrturbative in A, we will deter- 
mine a general expression for a field that is close to being 
Gaussian. The distribution of maxima and minima can 
be calculated from the joint distribution of h z , h zz and 
h zz * , which can again be determined given the cumulants 
of these variables. We will assume that the third order 
cumulants are of order A, while higher order cumulants 
are of higher order. (The formula is valid for the KPZ 
problem, since the assumption on the cumulants is sat- 
isfied for any function that evolves nonlinearly from an 
initially Gaussian signal, as long as the quadratic terms 
have coefficients of order A, and for times that are not 
too large.) The relevant cumulants up to third order are 
named in Table ITT1 

The expression for the imbalance An is derived by writ- 
ing the joint probability distribution in terms of the cu- 
mulants, as before. We then set h z — and integrate 
over the range of h zz and h zz * that defines the maxima 
and minima, to find 



An 




(7) 



Equation ^ is a general result. In order to apply it 
to the KPZ equation, we first change variables to 



u(r,t) = 



2v 



A 

h(r,t) 



exp 



h(f,t) 



4// 



h{r,tf 



(8) 



Note that u is a monotonic function of h, so u has the 
same profile of maxima and minima as h. This new field 
satisfies the heat equation whose general solution is 



where r = k^vt. The validity of this equation is illus- 
trated in Fig. |1J which shows an excellent agreement be- 
tween theory and numerics. The imbalance starts out at 
zero because the initial choice for h is Gaussian. On the 
other hand, after long times, this expression decays back 
to zero. The reason for the decay is that, at long times, 
it(r, t) involves an average over a larger and larger win- 
dow, so by the central limit theorem, it starts to acquire 
Gaussian statistics characterized by a vanishingly small 
imbalance between maxima and minima. 

For early times, we can make an expansion in t, valid 
for an arbitrary power spectrum, which gives 



An 



A I 6 



1 



V 9 V 7T K- 



--(2K 2 K 6 - 3Xf )(i/*) 2 + 0(t 3 ). 



(11) 

Note that for the Gaussian power spectrum, featured in 



Eq. (10), the second order term happens to vanish. 



The agreement for the entire range of times is special 
for the KPZ equation: for more general equations, the 
agreement would break down at sufficiently late times, 
since the nonlinearities eventually grow exponentially; 
the exact transformation of the KPZ equation to a linear 
equation implies that the nonlinearities remain bounded. 



Appendix A: Determination of g(z) 

The probability distribution g(z) can be derived using 
a method similar to the one outlined in [9 . Consider a 
fixed point r. We wish to know the probability density 
that at this point we have Hg = z (to avoid confusion 
with the derivatives of Hg, we shall write H from now 
on) given that it is a minimum. The conditions for this 
can be written in terms of derivatives of H, namely, H x = 
H y = 0, defines a critical point while H xx H yy — H xy > 
and H xx + H yy > distinguishes a local minimum from 
a saddle or maximum. First let us determine the joint 
distribution of these six variables (H and its derivatives) , 
which form a set of correlated Gaussian variables. The 
joint probability distribution - . - for any such 
set is completely determined by the correlations between 
the variables: 



u{f,t) a J d 2 r>G(r,r>,t)(Hr',0) + ^%',0) 2 ), (9) 

where G(f,r',t) denotes the Green's function. 

The correlations listed in Table 2 can now be deter- 
mined from the distribution of h(f, 0), leading to an ex- 
pression for An(i) that is valid over an arbitrary time 
span provided that A is small. Analytical results can 
be obtained for a few convenient choices of the power 
spectrum of h(r, 0). For example, if we take a Gaussian 
spectrum, A(k) 2 ~ exp(— k 2 /2kg), we find 

A» = V 16 ^ 1 + 4 ^ , (10) 
v V3^(l + 2r) 3 (l + 6r) 4 V ' 



p(6, ...,£„) = (2tt det Cy n ' 2 exp I - £ r /Ci, j , 

(Al) 

where = is the matrix of correlations between 

the variables. 

Correlations between H and its first and second deriva- 
tives can be expressed in terms of the first three mo- 
ments (K , K4) of its amplitude spectrum. By dif- 
ferentiating the Fourier expansion of H we find that 
(H 2 ) = (H 2 ) = ^K 2 , and likewise that the variances 
of the second derivatives are proportional to K4. The 
only variables among the six that are correlated to one 
another turn out to be H, H xx and H yy , with (HH XX ) = 
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(HHy 



-K 2 /2 and (H xx H y 



After re- above, the monstar fraction depends only on the joint 



trieving the probability distribution, we set H = z, 
H x = H y = and integrate over H xx , H yy and H xy 
(over the domain defining a minimum). The Jacobian 
determinant \H xx H yy — H xy \ must be added The 
probability density we have calculated so far reflects the 
chance that H x and H y are close to zero at the point 
r~o (there is a vanishing chance that they are exactly 0). 
However, we want to find the probability of the reverse 
situation, that H x = H y = exactly at a point within 
a small range of tq (since we are looking at the distri- 
bution of extrema in the plane). The ratio of the two 
probabilities is given by the Jacobian determinant. The 
final answer reads 




where A = 



- 2 (3-2A) er f c 



A(l - z 2 )e~2 z erfc 



: V'3A(l-A)ze" 5 ^T : 



A" 



A 



2(1-A)(3-2A) 



A 



2(1 -A) 



(A2) 



K K 4 



Appendix B: Complex notation and cumulants 



distribution of H zz , H zzz , H zzz * and their conjugates, 
and the properties of a Gaussian field are determined 
by covariances, there are only a few variables the mon- 
star fraction can depend on. These are a = (\H ZZ \ 2 ), 
t = {\H ZZZ \ 2 ), and r' = (\H ZZZ ,\ 2 ). The latter two are 
equal thanks to translational symmetry, and since there 
is no dimensionless function of a and t, olm must be a 
constant. 

Similar arguments can be used to show that the shift 
in olm for the non-Gaussian h = H+ef(H) depends only 
on (f"'(H)) to first order. First, we need to determine 
all the cumulants of h zz , h zzz and h zzz * up to first order 
in e. At this order, we will show that most of the cumu- 
lants vanish. We are faced with cumulants of the form 
C n (Dih, D2I1, . . . , D n h), where the operators Di repre- 
sent two or three derivatives. Each operator acts on h(r) 
at the same point r. For the moment, we will consider 
each operator Di to act on a different point f\. This al- 
lows us to bring all the operators outside the cumulant. 
This yields D\ . . . D n C n (hi, ■ ■ ■ , h n ), where hi — h(fl). 

Since we are only working up to first order, we can 
expand the cumulant as 

C n (hi, . . . , h n ) — C n (Hi, . . . , H n ) 

+ eC n (f(H 1 ),H 2 ,...,H n ) + ... 
+ eC n (H u H 2 ,...,f(H n )), 

(B2) 



Introducing the complex variables z = x + iy and z* 
allows one to use isotropy to calculate the probability 
distribution very efficiently. The isotropy causes many of 
the cumulants of the distribution to vanish - this explains 
why the monstar fraction is always the same for Gaussian 
fields when they are isotropic and why the shift of the 
monstar fraction depends on the nonlinearity f(H) in a 
simple way. 

The isotropy implies that a moment like (h zz h zz * z *) 
does not change when the field is rotated by an angle a. 
On the other hand this rotation transforms z —> ze la 
and z* — > z*e~ la , which in the given example would 
introduce an extra factor e la in the moment. Since the 
moment cannot change, it must be vanishing. In general, 
a moment can only be nonzero if the number of z and z* 
derivatives match. 

In these complex variables, the definition of an umbilic 
is given by h zz = 0. Monstars are distinguished from 
lemons and stars with the conditions 



\ 2 - \h zzz \ 2 > 0, (Bla) 



271/1, 



18\h 



?\h .I 2 

zzz I I ,l zzz I 



^{h , zzz h zz * z * + h z * z * z * h zzz * ) > 0, 



(Bib) 



which translates in complex notation the condition de- 
rived in Ref. [8] for an umbilic to be a monstar. 

First, we review why the monstar fraction is a constant 
for a Gaussian field H. Since by the conditions outlined 



which consists of one leading order term and n first order 
terms for which we can apply 



(f ( - n ^ 1 \H 1 ))(H 1 H 2 )(H 1 H 3 ) . . . (HA 



(B3) 



In this expression, the operators D 2 through D n can eas- 
ily be reinserted. For the operator Di the product rule 
needs to be applied, which in principle gives rise to a 
lot of terms. Remember though that, after setting all 
fl equal again, each moment can only be nonzero if the 
number of z and z* derivatives match. This criterion 
kills most of the terms. For example, if we consider 
the cumulant C^(h zzz * , h zzz * , h z * z *), we encounter the 
term d ltZZZ *{f" {^^(HxH^iHxH^^). In the prod- 
uct rule, after setting f[ equal to the other fl, the only 
nonzero term is (f"(H))(H z *H zzz ,)(H zz H z , z »). 

We can now also see that there are no cumulants C n 
with n > 4 which are nonzero up to first order in s. This 
is because when we apply the above recipe, we have n—1 
moments of the form (HiDiHi). The operator D\ has 
only three derivatives at most, therefore, in every term 
in the product rule there must be at least one moment 
of the form (HDiHi). However, the variables that we 
consider, h zz , h zzz , h zzz * and their conjugates, all do 
not have the same number of z and z* derivatives. That 
means that this moment must be zero, and hence the 
entire term too. 
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